function out = getImpulseResponsesCondStd(params,outEstim,IRFlength,shockSizeScalar,xfGIRF,xsGIRF,shockNum,plotGIRFsON)
numSim    = 250;
varToSave = {'$l_t$','$c_t$','$\pi_t$','$r_t$','$r^{(40)}_t$'};
%varToSave = {'$l_t$','$w_t$','$\pi_t$','$r_t$','$r^{(40)}_t$'};

%% Impulse response functions
reduceStatesON = 1;
selectY = 1:6;
%nn = sum(selectY~=0);
[model,errorMes]  = perturbationDSGEmodel(params,outEstim.model.orderApp,outEstim.setupFilter,selectY,reduceStatesON);
model.inclSurveys = 0;

% Test the size of xfGIRF and xsGIRF, the point around which the GIRFs are
% computed
nx  = size(outEstim.model.hx,1);
if ~isempty(xfGIRF)
    if size(xfGIRF,1) ~= nx
        error('xfGIRF must be of size nx*1')
    end
end
if ~isempty(xsGIRF)
    if size(xsGIRF,1) ~= nx
        error('xsGIRF must be of size nx*1')
    end
end
eta      = model.eta;
rng(1,'twister')
shocks          = randn(model.ne,numSim);

% Simulating the GIRFs
dimy            = model.ny+length(model.maturities);
IRFcondStd_sim  = zeros(dimy,IRFlength,numSim);
for s=1:numSim
    if mod(s,50) == 0
        disp(['Draws so far: ',num2str(s)])
    end
    xf_cup       = zeros(model.nx,IRFlength);
    xfShocks_cup = zeros(model.nx,IRFlength);
    xs_cup       = zeros(model.nx,IRFlength);
    xsShocks_cup = zeros(model.nx,IRFlength);
    xrd_cup      = zeros(model.nx,IRFlength);
    xrdShocks_cup= zeros(model.nx,IRFlength);
    xf_cu        = xfGIRF;
    xfShocks_cu  = xfGIRF;
    xs_cu        = xsGIRF;
    xsShocks_cu  = xsGIRF;
    xrd_cu       = zeros(model.nx,1);
    xrdShocks_cu = zeros(model.nx,1);

    
    for k=1:IRFlength
        % with shock
        shocksUsed = shocks(:,s);
        if k == 1
            shocksUsed(shockNum,1) = shockSizeScalar;
        end
        xfShocks_cup(:,k) = model.hx*xfShocks_cu + eta*shocksUsed;
        xsShocks_cup(:,k) = model.hx*xsShocks_cu + model.HHxxtil*kron(xfShocks_cu,xfShocks_cu) + 0.5*model.hss;
        xrdShocks_cup(:,k)= model.hx*xrdShocks_cu + 2*model.HHxxtil*kron(xfShocks_cu,xsShocks_cu) + ...
            model.HHxxxtil*kron(xfShocks_cu,kron(xfShocks_cu,xfShocks_cu)) + 3/6*model.hssx*xfShocks_cu + 1/6*model.hsss;
        
        % the benchmark without shock
        xf_cup(:,k) = model.hx*xf_cu + eta*shocks(:,s);
        xs_cup(:,k) = model.hx*xs_cu + model.HHxxtil*kron(xf_cu,xf_cu) + 0.5*model.hss;
        xrd_cup(:,k)= model.hx*xrd_cu + 2*model.HHxxtil*kron(xf_cu,xs_cu) + ...
            model.HHxxxtil*kron(xf_cu,kron(xf_cu,xf_cu)) + 3/6*model.hssx*xf_cu + 1/6*model.hsss;
        
        % Update of states
        xfShocks_cu = xfShocks_cup(:,k);
        xf_cu       = xf_cup(:,k);
        xsShocks_cu = xsShocks_cup(:,k);
        xs_cu       = xs_cup(:,k);
        xrdShocks_cu= xrdShocks_cup(:,k);
        xrd_cu      = xrd_cup(:,k);
    end
    % The controls
    pruningOn    = 1;
    %tmpShocks    = simConStdY_linear([xfShocks_cup;xsShocks_cup;xrdShocks_cup],model,pruningOn);
    %tmp          = simConStdY_linear([xf_cup      ;xs_cup      ;xrd_cup],model,pruningOn);
    tmpShocks    = conStdY_sim([xfShocks_cup;xsShocks_cup;xrdShocks_cup],model,pruningOn);
    tmp          = conStdY_sim([xf_cup      ;xs_cup      ;xrd_cup],model,pruningOn);
    
    % IRF in deviation from steady state value
    %IRFcondStd_sim(:,:,s) = (tmpShocks.conStdYDevSS(:,:)-tmp.conStdYDevSS(:,:));
    
    % IRF in levels
    IRFcondStd_sim(:,:,s) = (tmpShocks.conStdY(:,:)-tmp.conStdY(:,:));
end
IRFcondStd = mean(IRFcondStd_sim,3);

%% Computing GIRFs for the conditional means
% For the elements in g - i.e. without bond yields
shockSize      = zeros(model.ne,1);
shockSize(shockNum,1) = shockSizeScalar;
[IRFcondMeanMacro,IRFx] = GIRFPruning_v2(model.gx,model.gxx,model.gss,model.gxxx,model.gssx,model.gsss,...
    model.hx,model.hxx,model.hss,model.hxxx,model.hssx,model.hsss,model.eta,...
    shockSize,model.vectorMom3,xfGIRF,xsGIRF,model.orderApp,IRFlength);

numYields = size(model.byx,1);
IRFcondMeanYields = GIRFPruning_v2(model.byx,reshape(model.byxx,numYields,nx,nx),...
    model.byss,reshape(model.byxxx,numYields,nx,nx,nx),model.byssx,model.bysss,...
    model.hx,model.hxx,model.hss,model.hxxx,model.hssx,model.hsss,model.eta,...
    shockSize,model.vectorMom3,xfGIRF,xsGIRF,model.orderApp,IRFlength);
labelYields = cell(1,numYields);
for i=1:numYields
    labelYields{i} = ['$r^{(',num2str(i),')}_t$'];
end

IRFcondMean = [IRFcondMeanMacro;IRFcondMeanYields(model.maturities,:,:)];

labelyFull  = [model.labely labelYields(model.maturities)];
%% Selecting the variables we want to save
SelectVar = zeros(1,length(varToSave));
for i=1:length(varToSave)
    SelectVar(1,i) = find(strcmp(labelyFull,varToSave(i)));
end
% Updating the labels
labelyPlot = [labelyFull(SelectVar),labelyFull(SelectVar)];
labelTmp   = labelyFull(SelectVar);
for i=1:length(varToSave)
    tmp = char(labelTmp{i});
    labelyPlot{length(varToSave)+i} = ['$std_t(',tmp(2:end-3),'_{t+1})$'];
end

%% Plotting
if plotGIRFsON == 1
    nyAll    = size(IRFcondMean,1);
    scalingY = ones(nyAll,1)*100;
    IRFyPlot = zeros(nyAll*2,IRFlength,3);
    IRFyPlot(1:nyAll,:,model.orderApp)         = IRFcondMean(:,:,model.orderApp);
    IRFyPlot(nyAll+1:2*nyAll,:,model.orderApp) = IRFcondStd;
    % Scaling
    IRFyPlot(1:nyAll,:,:) = repmat(scalingY,1,size(IRFyPlot,2),size(IRFyPlot,3)).*IRFyPlot(1:nyAll,:,:);
    IRFyPlot(nyAll+1:nyAll*2,:,:) = repmat(scalingY,1,size(IRFyPlot,2),size(IRFyPlot,3)).*IRFyPlot(nyAll+1:nyAll*2,:,:);
        
    
    labelyFullwithStd = [labelyFull,labelyFull];
    labelTmp   = labelyFull;
    for i=1:length(labelyFull)
        tmp = char(labelTmp{i});
        labelyFullwithStd{length(labelyFull)+i} = ['$std_t(',tmp(2:end-3),'_{t+1})$'];
    end
    
    IRFxPlot = zeros(model.nx,IRFlength,3);
    shockSize             = zeros(model.ne,1);
    shockSize(shockNum,1) = shockSizeScalar;
    plotImpulseResponseCondStd(IRFyPlot,...
        IRFxPlot*100,labelyFullwithStd,model.labelx,shockSize,model.orderApp);
end

%% Saving output
out.IRFcondMean = IRFcondMean(SelectVar,:,model.orderApp);
out.IRFcondStd  = IRFcondStd(SelectVar,:);
out.labelyPlot  = labelyPlot;
out.shockNum    = shockNum;
out.numSim      = numSim;
